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Abstract. We present a new study of the evolution of the Carina dwarf galaxy that includes a simultaneous 
derivation of its orbit and star formation history. The structure of the galaxy is constrained through orbital 
parameters derived from the observed distance, proper motions, radial velocity and star formation history. The 
different orbits admitted by the large proper motion errors are investigated in relation to the tidal force exerted 
by an external potential representing the Milky Way (MW). Our analysis is performed with the aid of fully 
consistent N-body simulations that are able to follow the dynamics and the stellar evolution of the dwarf system 
in order to determine self-consistently the star formation history of Carina. We find a star formation history 
characterized by several bursts, partially matching the observational expectation. We find also compatible results 
between dynamical projected quantities and the observational constraints. The possibility of a past interaction 
between Carina and the Magellanic Clouds is also separately considered and deemed unlikely. 

Key words. Carina dwarf galaxy, Local Group, dwarf galaxies, chemical evolution 
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1. Introduction 

It has been e x tensively reviewed by many authors (e.g., 
iGrebell (fl999h : lvan den Berghl fll999h : llVlateol (Il998l )) that 
dwarf spheroidal galaxies (dSphs) are the most common 
type of galaxies in the Local Group (LG). Their origin 
and evolution are focal points because of the key role 
that they play in the context of the hierarchical growth 
of structures in the Universe, and because of the natural 
interest that these peculiar systems engender as exam- 
ples of dark matter dominated systems. In the context 
of dynamical research on interacting systems, dSphs also 
play a crucial role as stellar systems found in the central 
regions of galaxy clusters and groups, thus becoming po- 
tential subjects of strong gravitational interactions. This 
morp hology-positi on relation is also observed in our LG 
(e.g.. lGrebellll999l) . 

Tidal interactions may leave a variety of struc- 
tural signatu r es (se e, for instance 



Mufioz et al 



Peharrubia et al 



Read et al 
(I2008L 




In the LG, Sagittarius is the best-known example of a 
currently disrupting dwarf galaxy, but based on photo- 



metric studies tidal interactions have also been claimed 
for ot her dSph s such a s U rsa Minor, Sculptor and Carin a 



Eskridgd dl988alfbh: llrwin fc Hatzidimitrioul (Il995l) 
iKuhn et al.| jl996h: ISmith et all (Il997l): IMaiewski et al 



)5, 
)3) 



200 al2000h:lMartmez-Delgado et al.1 (l200ll):lPalma et al 



2003); Mufi oz et al 
Westfall et all <l2006h: 



(2005) 



Walchcr et al. 



Wilki nson et all (120 04); Sohnetal 



(2003) 



( 20071); iMunoz et all (|2006h : IMaiewski et all (|2003l ); 
Stri gari et all ( 20081 )). Draco on the other hand, appears 
to be an example of an un disturbed dSph despite it s 



proximity to t he M ilky Way (jOdenkirchen et all (|200ll ): 



Kles sen et al. I (l2003l )). Tidal interactions and their effects 



have also b een extensively explor e d in N-body simula - 
tions (e.g. 



Pasetto et al 



Johnston et all (jl999lk iMaver et all (|2002l ); 
20031) 1 In this context, a general bimodal 
density profile can be used to describe the systems with 
an inner stellar population and a shallowly decreasing 
density profile in the outer regions as a consequence of the 
tidal interaction. But while the external density profile is 
generally attributed to the tidal interactions, it is unclear 
whether the inner regions of dSphs are strongly tidally 
influenced, or whether a large M/L ratio can dampen the 
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tidal gr avitational shocks 



1997); Gomcz-Flcchoso et al 



see iBurkertl (Il997h 
(Il999l): 



Kroupa 2.2. Observations 



Klessen fc Zhao 



20021): iGomez-Flechoso fc Martmez-Delgadd (I200.C 2.2.1. Proper motions, distance and radial velocity 



Clevna et al.l 



Gilmore et al 




Walker et al.l (l2006aUbh : 
The radial velocity disper- 
sion profile is a powerful tool for the investigation of 
multi-component self-gravitating systems, since it is 
sensitive to the dark matter distribution and accessible to 
observations (e.g. Pasetto et al. (2010)). Recently radial 
velocity data became available to track the line-of-sight 
velocity dispersion as a function of radiu s for many LG 



Walker et all ^OOBaUbTl 



dwarf g alaxies (e.g. iTolstoy et all (120041): IWestfall et ail 
(120061); I Wilkinson et al.l fo004; iMuhoz et al.l (120051); 



iMunoz et al 



Sohn et al 
(120061 ) 



(200 71): iKoch et al 



Walker et al.l (120091^ 



f|2007c 

and these 

the dwarf galaxies kinem a tic status (e.g. iKroupa 



data permit a more detailed mod eling of 



dl997l): iKieyna et al.l (Il999l): iKazantzidis et all (|2004l ); 
Read et all (|2006l) : IPefiarrubia et al.l (|2008l.l2009h . Pasetto 
et al 2009). For the same approach to multi-component 
dSphs in the cosm ologi cal ACMD context see, e.g., 
Stoehr et all ([20021 ) and lHavashi et all d2003l) or. in a 



purely dynamical context see e.g. ICiotti fc Morganti 



ery 

( 20091) . Pasetto et al. (2010), where the primordial stellar 



population is embedded inside an extended dark matter 
halo. 

In our current study, we concentrate on the Carina 
dSph galaxy, for which an extensive body of observational 
data is available from the literature. We explore the orbital 
evolution of this dwarf. The paper is divided in two main 
sections: Sect. [2] reviews the observations on Carina and 
presents a preliminary orbit investigation in a point-mass 
approximation, while Sect. [3] presents the N-body simula- 
tions and comparisons with the observations. Conclusions 
are presented in Sect. |U 



2. The family of orbits for the Carina dwarf galaxy 

2.1. Initial condition vs final condition 

In studying the orbits of the Carina dwarf galaxy, aspects 
of its actual dynamical and kinematic condition, as well as 
the history of its stellar populations and chemical enrich- 
ment, must be taken into account. The incompleteness of 
available data, together with their uncertainty or errors, 
combined with the limits of our analysis tools, increase 
the difficulties we have to overcome in order to produce 
a realistic model for Carina's orbit. Here, the complete 
amount of available data for this galaxy is too large to be 
completely reviewed, thus only the more salient charac- 
teristics are reported for the sake of brevity. 



1 A slightly differ ent ap proach can be found in 

iPiatek fc Prvorl (|l995l ) and lOh et al.1 (|l995l ) where the 
authors span, in a specialized context, a wide range of eccen- 
tricities for exceedin gly low-mass dwarf gal axies without dark 
matte r (but see also Klessen fc Zhaol l|2002l ) and lKlessen et all 
l|2003l )). 



The present observational data on the Carina dSph 
make it necessary for us to discuss families of orbits 
instead of a true single orbit. Here, we devote our 
main efforts to reducing the number of plausible or- 
bits using all the available constraints known to date. 
We start with a determin a tion o f the proper motion 
in Table 1 of Metz et ail |2008) based on advanced 



charge transfer inefficiency (CTI) correction for the Space 
Telescope Imaging Spectrograph (STIS): (fi a cos<5, 11$) = 
(+22 ± 13, +24 ± 11) mas ■ yr^ 1 . We proceed with the re- 
duction t o an inertial galactocen t ric refe rence system So, 
updating IJohnson fc Soderbloml dl987l) to J2000 as al- 



ready done by e.g., IPasetto et al.l ( 20031 ). We extensively 



describe the orientation of the velocity space directions for 
this inertial reference system centered on the MW, So, in 
Appendix A. The resulting velocity vector for Carina is 
v C ar = { + 113 ±52, -14 ± 25, +44 ±56} km -s" 1 . Proper 
motions for Carina have also been recently derived by 
Walker et al. ( 20081 ) with error bars compatible with the 
values we adopted here. 

The first determination for the distance of Carina 
comes from the RRLyra e-b ased distance modulus 



estimation by ISaha et al 
cent values were 



and 



more re- 



19861), 

determined by ISmecker-Hane et al 



( 19941): iHurlev-Keller et all dl998l) : Mighelll dl997l) and 
McNamaral (|l995l ). 



Measurements of the radial velocit y were already mad e 



in an older work on carbon stars by ICook et al 



where the authors deduced a radial velocity dispersion 
of 6km • s _1 . The actually value we a ssumed for the ra- 



dial velocity is v r = 224 (±3) km • s _i (|Mateo et al.lll998 



Mould et al.lll990l) 



Methodology for the dwarf galaxy simulation: 

We simulate a dwarf galaxy composed of three compo- 
nents: a dark matter component that dominates the total 
mass, a gas component and a stellar component. We as- 
sume that the preexisting old stellar population has settled 
after an initial violent-relaxation phase of the dark halo 
which created the potential well within which the bary- 
onic component collapsed. The corresponding burst of star 
formation created the oldest stellar population that is ob- 
servable through the analysis of the observed colour mag- 
nitude diagrams (CMD) (see next section). We assume a 
time-scale of a few Gyr for the realization of this initial 
3-component dwarf galaxy model and its subsequent set- 
tling down in an orbit around the Milky Way (MW) with 
an initial orbital energy E™^. We start our simulations 
after this initial evolution has already taken place, i.e., we 
consider a dwarf galaxy in which the old population has al- 
ready formed and that is already in orbit around the MW. 
We now consider the subsequent orbital evolution cover- 
ing the last = 9 Gyr. Hence we assume that the observed 
proper motions represent are indirectly related this initial 
orbital energy E 1 ^. There is no reason to investigate the 



S. Pasetto et al.: Carina dwarf galaxy 



3 



orbit prior to this initial time, because the system has still 
not settled down in thermodynamic/virial equilibrium and 
the concept of a barycenter for a complex of blobs of gas 
and stars is still meaningless. Similarly, any old preexisting 
stellar population, by experiencing a collisionless violent 
relaxation process, loses memory of its initial phase-space 
distribution that, as a consequence, is completely unre- 
lated to the currently observable kinematic properties of 
Carina (see following sections for a review of other possible 
scenarios of formation for Carina) . 
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2.2.2. Milky Way model 

In our simulation we included an external galaxy model 
resembling that of the Milky Way. Nevertheless, consider- 
ing the large Galactocentric distance of the dwarf galaxy 
we are going to analyse and the large uncertainties in the 
MW models, we will model the MW with a few simplified 
assumptions: 

1. Milky Way potential. The Milky Way will be sim- 
plified as a vector field across which we will orbit first 
a point mass model of the Carina dwarf galaxy (in 
Sect. 12.31) . Then, once we have reduced the plausible 
orbits suitable to our aims, we will orbit a full self- 
consistent N-body model of the Carina dwarf galaxy 
considering gravity, star formation and feedback pro- 
cesses of this 'live'-satellite. We tested different poten- 
tials representing the MW in order to obtain results 
independent from our parametrization. Nevertheless, 
we want to stress that the large errors involved in the 
observational proper motion determination in general 
do not permit us to easily constrain the MW potential 
from the orbit determination of its satellites. The fully 
analytical parametrization is based on a logarithmic 
halo 



& h (R,z,t) = kflog R 



z 

" 2 + k 2 



32 1 kl + 



(1) 



for th e disk we can choos e a po tential formulation as 
(after iMivamoto fe Naeail (|l975h ) 



$d (R,z,t) 



and for the bulge 



GM d 



R 2 + (hi + v^Tfcj) 



(2) 



$ b (R,z,t) = - 



GM h 



s/R 2 + z 2 + k 6 



(3) 



where the dependence on the position is made explicit 
by the cylindrical parametrization (O, R, z) of So while 
the dependence on the time t is implicit in the pa- 
rameter ki = ki (t) for i = 1,...,6 and on the total 
masses M d = M d (t), M h = M h (t) for the disk and 
the bulge respectively. G is the gravitational constant 
G = 6.673 x 1Q- 11 m 3 kg- 1 s~ 2 . Particular attention 
is paid to the tuning of these parameters in order to 
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Fig. 1. Rotational curve of the growing halo of the MW 
potential. The shadow represent the zone of gained mass 
for accretion in the outer part of the halo. 



match the nowadays known constrains on the MW po- 
tential. The adopted parameter are listed in Table [TJ 
The evolution of the rotation curve for the radial range 
of Galacocentric distances of interest for the Carina's 
orbits is shown in Fig. [T] The upper blue thick line 
is e.g. compa t ible w ith the recent determination by 
Gnedin et al. (2010) even if larger margin of error in 



the rotation curve is a ctually permitted in reltaion to 
the Su n location 
( 20091) : 



e .g.. iReid et al.l (|2009t ): iBovv et al 



Sirko et al 



(2004^ 



dif- 



We stress that different studies have found 
ferences for the MW mass interior to the zone 

of releva nce for our orbit integration. E.g., 

Kochanekl (| 19961) obtained M tot (r < lOO kpc) = 
7.5 x lO n M , and lDehnen fc Binnevl (jl998l) obtained 
M tot (r < lOOkpc) = 7 x lO n M . In addition, many 
different appr oache s can produce mass estimates for 
the MW, e.g. IZhaol (|l998l ) used the velo ci ty dis tribu- 



tion of the MW satellites, iRobin et al.l ([20031 ) used 



the escape ve l ocity from the Local Standard of Rest, 
and Combesl ( 19991 ) fro m timing argum e nts. S imilar 



results were found by ISakamoto et al 



the kinematics of Galactic satellites, 



(|2003l ) from 
halo g l obular 



clusters and horizontal branch stars. iLin et al.l (|1995[ ) 
obtained M tot (r < lOOkpc) = 5 x 10 n M^ from the 



dyna mics of the Magellanic Clouds, IVallenari et al 
(|2004l) found Mtot (r < lOOkpc) = 8.8 x 1O U M from 
proper motion modelling of stars in two deep fields 
toward the North Galactic Pole. Mass estimate based 
on the moments of the Jeans equation are often 
sensitive to t he an isotropy profile of the MW (e.g. 
Brown et al.l (l2010l)'l leadi n g to u nsecur ed results, e.g 



(2006) 



compa re iBattagli a et al 



or 



Wat kins et al 



Dehnen et al 



2005) and 

(|2010l) that which quoted 
a value for M tot (r <_300kpc) « 1.4 x 10 12 M Q and 
recently IXue et al.l (|2008l) found M tot (r < 60kpc) 
4.0 x lO n M . 

Milky Way accretion rate. In these last 9 Gyr of 
evolution, we expect that not only the dwarf galaxy 
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Table 1. Parameters adopted as starting and final values for analytical MW potential of Eqn. (JTJ) , @ and ([3]). 



^ini 


M b 
[M ] 


M d 

[M ] 


fci 

\kms~ 1 ] 


k 2 
[kpc] 


k 3 


ki 
[kpc] 


fes 
[kpc] 


kc, 
[kpc] 


Bulge 

Disk 

Halo 


3.4 x 10 iu 


0.9 x 10 11 


120.8 


10 


1 


6.5 


0.26 


0.7 


tend 


Bulge 

Disk 

Halo 


3.4 x 10 1U 


1.0 x 10 11 


130.8 


12 


0.8 


6.5 


0.26 


0.7 



changed its physical properties, but that also the mass 
of the MW grew due to the continuous inflow of mat- 
ter. For simplicity we assume that 

— The MW has not experienced any major merger in 
the last 9 Gyr, 

— star formation processes are of minimal importance 
in the growth of MW. 

These simplifications are realistic for o ur time range 



of interest of roug hly = 9 Gyr (e.g., iGuo fc Whitd 
(2008;)). What is most interesting for us is the evolu- 
tion of the density distribution of the outer part of the 
halo, which is typic ally expected to be natter than the 
inner part (see, e.g, iHelmil (|2008l ) for a general review 



on this dichotomy) , because the outer halo is where the 
orbit of a dwarf galaxy is expected to pass through and 
evolve. We represent the external mass distribution of 
the MW with an external potential that is time depen- 
dent. 

All these determinations lead to a wide range of val- 
ues within which we want to choose our parameters 
when a growth factor is added in a time-dependent 
potential for the halo component. For the MW we 
assume that in the time range of interest, roughly 
the last 9 Gyr, the growth rate is mainly due to the 
minor mergers that are expected to be dominant in 
general for re dshift z < 1 a nd fo r galaxies as mas- 
sive as MW ijGuo fc White! (|2008h . Taking into ac- 
count the uncertainties in the final target amount of 
dark matter, and the m ass to light ratio of the MW 



as in 



Flvnn et al. (2006). a simple linear growth fac- 



tor of 10 - 25% is tested for the MW in the last 9 



Gyr o f evolution (z < 1.46 in redshift) (e.g. ICombes 
( 19991) : IGuo fc Whitel (120081) ). Higher values for this 
parameter can be found in literature, although the 
relevence of the MW mass growing factor is limited 
by the large uncertentaintly in the total mass of the 
MW (M tot {r < lOOfcpc)). This uncertainties has un- 
avoidable bad influence on our ability to determine the 
orbits of all the dwarf galaxies in the MW outskirts, we 
also test and take into account a gradual flattening of 
the inner halo, although there exist no data on Carina 
dwarf galaxy that are sensitive to this parameter : e.g. it 



can be proved, as already done in lPasetto et al.l (|200 



that the presence of the bulge is not relevant for the 
integration of the orbit of Carina (similar as in the case 



of Sculptor in iPasetto et al.l (|2003l) and the following 
sectio ns). The on l y obse r vational data, e.g. from the 
SDSS lSirko etaD (|2004al) : I Juric et al.l (|2008l ). refers to 
the first 25 kpc of Galactocentric distance, that are ex- 
tremly unlikely in the orbit of Carina (see Fig. [5] and 
H]in the following section). 
3. Milky Way triaxiality. The Milky Way poten- 
tial and the possibility if its triax iality are in - 
vestigated in different paper s, e.g., IHelmil (j2004 . 



Martmez-Delgado et al.l (|2007f) . Unfortunately the ap 



proach used in these studies is based on the pos- 
sibility to t rack the orbits with the help of tida l 
streams, e.g.lCole et al.l (|2008l ): lEvre fc Binnevl |2009l ): 
Koposov et al.l (|20*i*0 ) and thus is limited to a few ob- 



jects, mainly globular clusters or to the Sagittarius 
dwarf galaxy. Our target, the Carina dwarf galaxy, 
does not show long tidal streams that would permit us 
to constrain the MW potential triaxiality. Hence, de- 
spite the simplicity of implementing a tiaxialhalo MW, 
we will not introduce here such a potential since its in- 
fluence cannot be directly verified within the errors 
of our present work. Moreover, triaxial halo presents 
challenges to the stabil ity of the disks as shown in 
Kazantzidis et al. ( 2010h . the origin and evolution of 



the dark matter axis ratios c/a and c/b are still un- 
clear, and the relative orientations of the three axes 
are still a matter of debate. Nevertheless, we stress 
that the role of triaxiality is invoked as fundamental 
tool to investigate the orbi ts of the MW sate llites close 
to the galactic centre (e.g. lLaw et ah I (12009^ . 



2.2.3. Star formation history 



In the literature, many 
star formation history 
since its discovery by_ 



pioneering work of Mould fc Aaronso: 
studies demostrated the presence 



authors have investigated the 
(SFH) of t h e Car ina dSph 
Cannon et al.l _ Jl977l) and the 



19831) . These 
of mu lt iple popula - 



tions with different ages ( LAzzopar di et al.l ( 1985 



(Il997l): iHurlev-Keller etal 



Smecker-Hane et al.l dl994l): iMould et all (ll99(f~lMighel]| 



1986) 



(Il998l ): iMighel (Il990ah ). 
i pu blished a deep CMP for the cen- 
and iRizzi et al.l (2003) derived the 



Monelli et al.l (|2003aD 
tral zones of Carina, 
SFH from a more spatially extended sample based on the 
analysis of RGB stars. We can then assume that the star 
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formation in Carina consisted of an initial burst of star for- 
mation that formed the oldest stellar population, a second 
important burst around ~ 6 Gyr ago, and a more recent 
event around ~ 2 Gyr ago. These episodes of star forma- 
tion are clearly visible in Carina's CMD, and the existence 
of a fourth, even more recent episode has been proposed. 
What makes Carina special is that these episodes are rec- 
ognizable as distinct events with distinct main sequence 
turn-offs. 

2.2.4. Chemical evolution 



Constraints 
only 



on the Carina metallicity distribution 
recently from the sp ectroscopic studies 



Koch et al.l (2006): 



Koch et all (|2007al ): lAbia et al.l 3 



came 

or chemical evo l ution models 
lLanfranchi et al' ( 2006 ): 
(|2008l) : iKoch et all (|2008h V They generally show a quite 
spread in the metallicity distribution. We will consider 
the chemical evolution in more details in a separate work 
(Pasetto et al. in preparation). 

2.3. Point mass integrations 

Since ou r work i s bas e d on two full N-bo dy integra- 
tors (e.g. lBerczikl (|l999( ): ICarraro et al.l (|l998l) ) which can 
track simultaneously many different physical processes 
and produce results that are compatible with observa- 
tional constraints, one might reasonably question the util- 
ity and validity of point-mass approximations. We con- 
tend that there are multiple advantages in the use of the 
point-mass approximation as a preliminary study for orbit 
determination: 

1. The principal effect that can influence the orbit deter- 
mination is the mass loss that the dwarf suffers dur- 
ing its pericentric passages. This mass change cannot 
easily be parametrized because it depends on the in- 
tensity of tidal effects as well as on the initial internal 
structure of the dwarf. Thus we will use phase-space 
coordinates and velocities obtained from point-mass 
calculations to get the initial values of the Carina or- 
bit for the full N-body simulations only for the first 
1 Gyr of backward evolution, i.e. before the time of 
the first pericentric passage t v \. These initial guess 
values for the phase-space location of the barycen- 
ter of the N-body system span a full range of po- 
sitions in configuration space from apocenter (where 
we will make an approximate match with the obser- 
vations) to first pericenter (at the time t = t p i) and 
lead to different star formation histories even if they 
correspond perfectly to the same initial orbital energy, 
E olh (i) = = E oxh (t )Vt G [t ,t pl ]. As an ex- 

ample of the orbit indetermination resulting from the 
large error bars on the proper motion, we plot in Fig.[2J 
[3] and HI the orbit integrations, backward in the time t, 
for point-mass particles where initial values have been 
randomly sampled with uniform distribution inside the 
error bars for the velocities. 



The symmetries of the orbits remain evident. The ex- 
pected appearance of the orbit of Carina from the 
point-mass integration is roughly a rosette-like orbit 
with several pericenter passages. This observation per- 
mits us to make some simplifying assumptions. The old 
stellar disk population of the MW is at least as old as 
10 Gyr (e.g. lFreeman fc Bland-Hawthorn! (|2002l )). thus 
we can assume that the gravitational potential of the 
disk and the bulge of MW has not changed its gravita- 
tional influence on the orbit of Carina. This allows to 
us to limit the growth rate in the MW mass only to the 
halo component, keeping fixed the amount of mass in 
the disk and the bulge. A growing disk has nevertheless 
tested as well in our simulation as explained in Table 

m 

The rosette orbit symmetry reflects a time symme- 
try that permits us to backward-integrate the N-body 
gravitational system by simply changing the sign of 
the velocities. Nevertheless, the present structure to 
which we apply the orbital phase space and the ex- 
ternal gravitational influence of the MW, describes a 
possible initial dynamical structure of the system as it 
was after the oldest stellar population had formed. It 
is often not considered in dynamical studies, in these 
kind of works, that the stellar evolution history does 
not have this time symmetry that is assumed in the 
pure dynamical approach. Clearly, the gas consump- 
tion and the chemical evolution are unidirectional also 
for the completely collisionless treatment of the dy- 
namics. Thus the point-mass integration provides us 
with an initial volume of the phase space values to as- 
sign to a full N-body simulation model at the starting 
time, say t = to = 9 Gyr ago. 

Even though it is clear that the orbital plane deter- 
mination of full N-body simulation and a point-mass in- 
tegration in a spherical potential is supposed to be not 
too different between, we would like to show that this is 
also the case for the pericenter passages in this type of 
orbit. This is done in Fig. [SJ Here, we plotted the results 
of the point-mass integration for an orbit with a close 
pericenter passage for the set allowed by Carina's proper 
motion error bars. The same initial conditions were then 
evolved in a fully consistent N-body simulation. The re- 
sults were obt ained with two different codes: the first by 
iBerczik ( 19991 ) in a serial machine w ith special hardware 
dedicated (see [Spurzem et al.l (|2009l) ) and the second by 



ICarraro et all (|l998l ) on a parallel machine at the Juelich 
Super-computer Center. Both codes yeld comparable re- 
sults across the entire range of initial values. 

So why does our approximation work so well? If we 
consider the two galaxies MW and Carina, we have to 
treat them as extended bodies. In this case the accelera- 
tion of the center of mass of Carina is given by 



,Car' 



dt 2 



Car (x')V x ,<W (x') 



Pg 



where the integral is extended to the entire galaxy and 
we can no longer assume Vx^a/w (x') = const to take it 
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Fig. 2. Error propagation for the point-mass approxima- 
tion in the orbit of the Carina dwarf galaxy. We see the 
convolution surface for 10 5 orbits. Although the uncer- 
tainty of the orbits does not permit us to go back in time 
for more than a few Gyrs, it is evident that a restricted set 
of values can be assigned to the orbital energy from the 
path of the galaxy in its first Gyr of backward evolution if 
we require to lie inside this manifold. This permits us to 
provide initial values for the full phase space description 
of the orbits that we want to integrate in a full N-body 
simulation, starting from the apocenter or the pericenter 
position. The use of these two different starting points as 
initial position of the barycenter of the gravitational sys- 
tem, although it refers to the same orbital energy, leads 
to different star formation histories. 



-100 
y [kpc] o 






z [kpc] 



-100 



Fig. 3. Same as Fig. [5] but with an edge-on view. The 
box-type shape of the orbit and the permitted regions are 
evidenced. 



outside the integral sign. This means that extended bodies 
orbit generally with the same velocities as a point parti- 
cle if and only if V x ' ( l > (x') = const, i.e. if V X '<I>(x') is 
roughly constant within the body volume of the Carina 
dwarf galaxy and it can be extracted off the integral sign 



& 60 
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10.25 
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Fig. 4. For the previous orbit integration shown in Figures 
[5] and SI we computed the radial distance of the dwarf 
galaxy barycentre in the respect to the inertial reference 
system centered on the MW. Here the probability of the 
location of the dwarf over the look-back time has been 
computed. Larger white zones (i.e., high probability) are 
closer to the present time where the indetermination has 
a minor impact. The indetermination affects mainly the 
pericenter passages where we have few or no constraints 
to satisfy due to the dominating effect of the MW growing 
halo potential. 
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Fig. 5. Radial distance evolution test. Here we compare 
the N-body simulation (black dots) and the point mass 
integration (red solid line). The dashed line represent the 
integration without the bulge commented in Sect. 12.51 The 
lines at the bottom represent the difference between the 
red-line and the dashed line, while the solid black line 
represent the difference between the red line and the orbit 
computed considering the dinamical friction (see text for 
details and Sect. I2.5|) . 
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to obtain the point mass equation of motion. In general, 
this approximation cannot be applied to galaxies because 
they are neither bodies that could be considered small 
enough such e.g., in the case of a globular cluster or- 
biting around the central zones of the bulge, nor can 
they move in a sufficiently homogeneous gravitational field 
e.g., they can reach the inner zones of the MW galaxy 
where the bulge influence become dominant (see Sect. 
12.511 . Using the consideration on the MW potential given 
in the previous section we have taken care to check that 
the orbital range of Carina in the external potential of 
the MW satisfies the required approximation in every or- 
bit computation performed with the point-mass approxi- 
mation and in minimizing the Action (see next section). 
Thus we safely proceed further with the approximation 
mcar (t) = mcar = const and the point-mass approxi- 
mation in the next section where we develop a dedicated 
method to select a preferential initial condition for the full 
N-body simulations. In the N-body simulations explained 
in the next section mass loss effects will be self-consistently 
taken in consideration. 



2.4. Minimizing the action 

After these preliminary considerations it became evident 
that, although the the errors in the proper motions are 
quite large, at least the orbital plane can be inferred rel- 
atively easily from the orbit analysis. 

For any chosen value within the range permitted by 
the error bars of the proper motions, the orbital plane re- 
mains roughly determined. The range of orbits allowed by 
the error bar permits us to deduce that the orbits are more 
likely box-type. That results not from the orbit determi- 
nation which is integrated for too short a time range, but 
from the consideration of full sample of orbits integrated. 

The true orbit, is probably contained in this plane but 
we want to proceed by constraining the possible orbits fur- 
ther with an alternative approach. The approach of using 
a point mass for the orbiting dwarf galaxy can bias the de- 
termination of the pericentre and apocentre of the orbit 
evolution if the dwarf passes very close to the inner zones 
of the galaxy. Nevertheless, we do not expect it to influ- 
ence the plane determination of the orbit, since the Carina 
dwarf galaxy is orbiting mostly in the external regions of 
a growing MW halo, nor do we expect to it to affect the 
error estimates of the pericentre position to be so large 
as to invalidate the following considerations. To support 
these approximations we remember that when Carina is 
more massive (before the first pericenter passage) the MW 
halo still has its spherical shape, which is gradually con- 
verted into an ellipsoidal shape at t en( i with parameter 
changing in Table [1]) from 1 to 0.8. Based on Fig. @] we 
can prove that the probability that Carina intersects the 
disk of the MW is null and the probability of a pericenter 
lower than 60 kpc is only a few percent! 



The generic Lagrangian of the dynamical MW-Carina 
system can be written down as 



||*(t)|| -m C ar$Gal(x,t) 



(4) 



The Lagrangian is a function that lives in the tangent bun- 
dle of the configuration manifold Q of Fig.[2]or|3l thus it is 
defined in L : TQxR — > E (where K is the infinite time line 
parametrized by the time t) and its integral can be easily 
solved numerically. Thus, we ask ourselves: what is the or- 
bital plane within the error bars of the proper motion that 
minimizes the action? The principle of minimal action will 
lead us to possible stationary points for the integral of the 

Lagrangian. Hence, we proceed by minimizing the action 
o 

integral S = 5 f Ldt, subject to the constraints on the ini- 

to 

tial value that we deduce from the proper motion and the 
star formation history. We can, in fact, suppose that the 
episodes of star formation have been triggered by the tidal 
influence of the MW potential in the pericentre passages 
(the proof of this hypothesis will come only in the next 
section). Although in TQ the equations of motion are a 
set of first-order differential equations (i.e., through each 
point (x, x) of TQ passes just one solution of the equa- 
tions of motion) thus permitting us to easily integrate the 
action S, before performing a minimization of the action 
we can already restrict the space T a , Car ( )Q, i.e. the fibers 
(in the geometrical sense) of the Carina configuration po- 
sition at the initial time a:car(0). 

In terms of Newtonian dynamics, we already 
gain a more restrictive velocity space than V 3 = 
{v\v E M. 3 /\ v — cr v ^ v ^ v + cr v } if we exclude from 
our velocity space all the orbits that do not lead to a 
pericentre passage at the times t%, t% corresponding to the 
two bursts of star formation at t 2 ~ 6 Gyr and t\ ~ 2 
Gyr. Moreover, we do not pretend to know the precise 
moment for the time of closest appro ach, but we can pro- 



1,2 



e.g. 



Hurley- Keller et al 



vide t he time range Sti, 

(jl998l )L In general we expect 8t 2 > 5ti and both not 
to be smaller that 1 Gyr. The desired conditions we 
want to implement in our minimization tool can easily be 
achieved by requiring a minimum in the distance function 
d (t, v )\ Cax : R x V 3 M in the range [t x - ^,t x + 
and in the range [t 2 — %■ , t 2 



61; 
2 



d t d(t,v ) 



OA d t ,td(t,v o )\ t=i >0 



(5) 



Here clearly the continuity of the partial derivative with 
respect to the time dt is assumed. We underline immedi- 
ately here that, as evidenced by Eqn. ([S]), this condition 
does not require at all to exclude other pericentre passages. 
In this case of the manifold at T XCar Q at t = results in 
the plot shown in Fig. |5]as an initial velocity space. 

At this stage, here we assume here that the tidal field 
acts as a trigger of the star formation. The proof of this as- 
sumption will be provide later with a selfconsistcnt gener- 
ation of the star formation history of Carina in the context 
of the smooth hydrodynamic particles approach integrated 
in a full N-body simulation. 
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Fig. 6. T XCar (Q) at t = 0. The section provides the zones 
within which the condition expressed by Eqn. ([5]) is sat- 
isfied. This condition clearly represents a good constraint 
when compared with the whole space where we want to 
find the minimum of the action. The advantage is in the 
more restricted number of points that we have to sample, 
for which we have to integrate the differential equation 
and to solve the action integral as well as in the reduced 
computational time. The volume of the velocity space to 
sample in this figure is more than 1/5 smaller than the 
overall space without the constraints imposed from the 
star formation history. 



Nevertheless, by excluding those orbits that do not 
show any pericentre passage in the distance range required 
by the matching of Carina's SFH will reduce further the 
space where we have to minimize the action (and the com- 
putational time required for the minimiz ation!). We adopt 



the same methods of minimization as Pasetto & Chiosi 



( 20071) . We point out that the method applied here does 
not require any previous or bit parameteri zation as in 
a cosmological context (e.g., Peebles! ( 1989h V Moreover 
the alternative approach that comes from the Fourier se- 
ries decomposition used in st ellar spectral dynamics (e.g. 
Carpintero fc Aguilar ( 19981 )) cannot be exploited here 
because of the small number of orbital periods during 
which we can integrate the dwarf galaxy (no more than 
9 Gyr and an orbital period between 2 and 4 Gyr) could 
lead to wrong frequency determinations. Nevertheless, for 
this small number of periods the full numerical integra- 
tion of the action works fine. The only limitation of our 
method is the treatment of the mass loss in Eqn. (TJ| : an 
approach for a dwarf galaxy experiencing a strong dwarf- 
bulge interaction as e.g., in the case of Sagittarius, is not 
able to reproduce correctly the orbits we obtain with a 



full N-body simulation. This limits our methodology to 
the more distant dwarfs such as the Carina dwarf galaxy. 

With this approach, we can generate a the family of or- 
bits with the following values (obtained through the mini- 
mization process of the action for the Carina dwarf galaxy 
starting at tit, — 9.645Gyr): 



a=Car = (-73, -64, 3) kpc, 
UCar = (65, -115, 75) km -s- 1 . 



(6) 



These are the initial phase-space values for most probable 
volume of T XCai ( t[b=9 .645Gyr)Q for the galaxy at the initial 
time of our N-body simulation. 

After a full N-body simulation run, the previous initial 
conditions lead to the following present-day best fit values 
(values for the barycentre of the galaxy in a full N-body 
integration at tn, = 0): 



zcar = (22, 89, -34) kpc 

v C ar = (-136, 13, -45) km • s~\ 



(7) 



(values in the reference system So) which are fully compat- 
ible with the proper motion, radial velocity and distances 
that are actually observed. 

2.5. Robustness of the point mass solutions 

The robustness of the the point mass integration is 
the subject of many cosmological or pure dynami- 
cal studies but often without a self-consistent consid- 
eration of the star format io n processes (e.g ., recently 
Penarrubia fc Benson! (j2005l ); iLux et all (|2010l )). The re- 
lation between star formation and orbit determination 
is first time fully exploited here as key ingredient to 
constraints the Carina orbit determination. Nevertheless, 
once we have determined a family of solutions as in Eqn. 
(J6j> , and before generating a complete N-body integration 
that self-consistently considers the structural parameters 
of Carina and its star formation history, we check this fam- 
ily of orbits against the po ssible role that the d ynamical 
friction could have played ()Chandrasekharl [ l9 43) and the 
parameters of the potential. 



2.5.1. The dynamical friction 

In the last century, the perturbative methods applied to 
dynamical friction studies have achieved several important 
results in order to overcome the limitation of the original 
formulation in the celebrated article by Chandrasekhar 
(i.e. ma inly the isotropy of the s urrounding stell ar 
field, see iTremaine fc Weinberg! ((l98l : IWeinbergl (Il986l) ). 
Closely related to t his ap p roach is also the work o f 
Bekcnstein & Maozl (|l992t ): iNelson fc Tremaind (|l999h 
that related the dynamical friction to the fluctuation- 
dissipation theorem. We will follow the treatme nt pre- 
sented in the paper of lColpi fc Pallavicin i (1998) which 
contains a s spe cial limit the Chandrasekhar work (e.g. 
Kandrunl (Il983h ). In this formulation of the perturba- 
tive theory (extension of previous works of e.g., iKandrup 
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( 198ll ): ISeguin fc Dupraz] l| 1994 1 19961 )) the authors recol- 
lect and extend several important results in order to over- 
come the limitation of the Chandrasekhar formula allow- 
ing one to follow analytically the gravitational wake influ- 
ence, the tidal deformation and the shift of the barycentre 
of the primary galaxy. Moreover, this formulation unifies 
the local and large-scale interpretation of dynamical fric- 
tion and provides a treatment even superior than the full 
N-body simulations often l imited by the number of parti- 
cles (see lCofpi et"aH (|l999h ). What is more interesting for 
us, is that in this last paper the authors already proved 
that the role of the dynamical friction for the MW dwarf 
satellite is marginal, especially with reference to the more 
distant dwarfs, and in conflict with what would be ex- 
pected from the previous straightforward application of 
the Chandrasekhar formula. Finally, if we want to perform 
a simple check, under the approximation presented in the 
paper of lColpi et al. I (|l999h with the further simplification 
that comes from the constraint on previous investigated 
relevant pericenter distance at which the orbit of Carina 
is expected to evolve have its pericenter distance, we can 
write rewrite the equation of motion in a system, say Si , 
as 



rf 2 x (t) 
' dt 2 



-GNm-^-z [ d^p(r') + F A , (8) 

||x(t)|| Jr'<x(t) 



where G = 4.4926 x lO" 6 kpc' j Gry" 2 M0 1 is the grav- 

MMwMpar reduced maS g 



itational constant, 



/' 



M Car +M M w 

of Carina and MW, p the halo density profile and 
F a the back r eactio n force as presented in Appendix 
of IColpi et al. (Il999h . An alternative faster integration, 
which given the distance of Carina is in our case abso- 
lutely equivalent in its results, is based on a time de- 
pendent Coulumb logarithm In A ~ I n j ra '" ~ ii ^ j with 

e ~ Csee lWhitd (| 19761) : I Weinberg! (119891) ) used in the 
the infinite, homogeneous, non-self-gravitating collision- 
less background formulation by Chandrasekhar: 



F A = -4tt(GM) pin A Erf (a;) 



2.i 



V 

rvr 



where x = and V si the velocity vector of Carina. 

The application of this formalism to our family of or- 
bits for Carina confirms that the role of dynamical friction 
on our orbit determination is irrelevant as shown in Fig. [5] 
The orbit computed with Eqn. ([5]) almost precisely over- 
lap the orbit determined without the dynamical friction 
(Sect. 12.41) thus we plotted a the bottom of the distance 
between the two orbits with a thin solid black line. We 
stress for comparison, that the error in the knowledge of 
the present location of the Cari na dwarf galaxy is of about 
10 kpc in radial distance (e.g. Mateo et al.l (1998)). 



2.5.2. The MW potential 

The role of the MW potential is fundamental in evaluating 
the orbit of Carina. Recently lLaw et aT ( 20091) proved the 



necessity to fully exploit the triaxiality of the MW poten- 
tial in order to explain the tidal streams location and the 
orbit of Sagittarius. Unfortunately the situation for Carina 
is much more complicated. The error in its present location 
and the absence of long tidal streams that can somehow 
track the orbits (or in the contrary can constraint the MW 
potential) prompted us to apply a different technique tak- 
ing the star formation history, as a new constraint on the 
orbit. I.e., we are explicitly assuming that Carina's un- 
usual star formation history was influenced by encounters 
with the MW. 

As an illustrative example we point out here how for 
the Carina orbit integration the existence whole presence 
of the bulge (~ 1O 1O M0) can be neglected! Carina is orbit- 
ing in zones where the relative influence of the potential is 
mainly dominated by the halo and partially by the disk: 
^buigc (x, i)/$tot (x, t) < 0.1 for every point x and instant 
t of interest in our orbit determination. By using the best 
fit orbit determination parameters (see Sect. 12.41 Eqn. ([5])) 
in Fig. [5] we plotted with the dotted blue line the orbit 
distance computed without the bulge. As evidence in the 
bottom of the plot, the distance between this orbit and 
the red orbit for our best solution retains a value compa- 
rable with the error bar as large as the present day Carina 
distance uncertainty (10 kpc). This result shows that the 
role of the bulge is marginal in the orbit of Carina. 



3. N-body integration 

The search for the best initial conditions of the orbital 
parameters of Carina must include also the matching of 
the internal properties of the dwarf. To accomplish this 
task we performed a large number of simulations in the 
space of the the restricted initial orbital condition that 
surround the best fit values obtained from the minimal 
action in the previous section. 

In order to perform this analysis in an automatic way, 
we wrote a code that computes all of the desired properties 
(derived star formation history, chemistry and kinematic 
properties of the orbiting evolved object) centred on the 
non-inertial barycentre of the particles representing the 
stellar component of the dwarf, say 5*2. The barycentre 
of the moving dwarf's stellar component also has to be 
automatically computed to speed up the analysis of the 
dwarf's orbits. This is performed with an original method 
based on the down-hill simplex algorithm (explained in 
Appendix IB"]) . 

A detailed description on the implementation adopted 
for the star formation recipes, chemical evolution 
and dyn amics fo r this GRAPE-treeSPH cod e can b e 

" (l2003al lbl. 120021) : 
and references 



found in [B erczik 
Berczik fc Petrov 



(1999) 



20011) 



Berczik et al 
Berczik (2000) 



therein. Special purpose hardware (jSpurzem et al. (2009)) 
was used in order to carry out a large number of high- 
resolution chemo-dynamical simulations. 
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3.1. Initial conditions 

The structural properties were adop ted from an ex t ension 
of the isolated models presented in iPasetto et al. (|2010l) 



profile by th e star formation processes as evidence in 
Pasetto et al. ( 2010f ). Flatter DM density profiles have 



in order to include a larger set of scale length initial condi- 
tions (Pasetto et al. 2010 in preparation). Those isolated 
models then are evolved in the orbits around the MW se- 
lected as explained in the previous section. Star formation 
history, chemistry and dynamical properties of the Carina 
dwarf galaxy are self- consistently determined along the 9 
Gyr of evolution. 

The staring initial condition can be partially 
tuned with the help of the literature results for 
Carina. Color magnitude diagrams were studied in 



dHurlev-Keller et al.l 


19981: iMighelll Il997t 


1 

1992, 


1990b.a; 


Monelli et al.l 12006b 


a, 2004b!lalld. 2003a 


Ibl: iRizzi et al.l 


2004 20031) ). Spectroscopic measurements yielded chemi- 



cal abunda nces, e.g.. iKoch et all (|2006l 120081 ). the M/L 
ratio, e.g., IGilmore et al. (120071) and the dark matter 



distrib ution, e.g., IHavashi et al.l (|2003l ): IGilmore et al 
' 20071). From the analysis of the CMD of Carina (e.g. 



Hurlev-Keller et al.1 (Il998l) l we inferred that the oldest 
stellar population contributes no m ore than 15% of the 
stellar mass and from the work of Mateo et al ] (Il993l) 



we got a mass to light ratio in the V band of (x)y 
23 (xf ) ■ Thus from Table 3 of IHavashi et"aH (|2003l) 
we got L — 0.4 • 10 6 L Q , which leads to a stellar mass 
estimate of M « 23 • 0.4 • 10 6 M Q = 9.8 • 1O 6 M . If 
we now assume that the mass lost by tid al interactions 



could reach an order of 95%, following e.g.. IHavashi et al, 
( 20031 ) we get M Ca r = 2 • 10 8 M Q and with an intial dark 
matter contents at l east 5 times the baryonic fraction 
Hurlev-Keller et al.1 (|1998f ). The baryonic material con- 



sists of 15% stars, 85% gas. We thus adopted a starting 



value of M, 



Mi 



15 



gas - = 33 ' 106M © and M star = 

-,. — 5.9 • lO 5 Af for our reference model. These 
numbers can be considered as starting values for an orbit- 
ing dwarf galaxy, aroung which we want to find our best fit 
initial values. The final best matching initial condition are 
reassumed in Tabled] The models are represented by a 3- 
component self-gravitating system where the family of the 
mult i- gamma model s was adopted for the de nsity profiles 
(e.g. iDehnenl ()l993l ): iTremaine etaH (|l994l) ). The scale 
parameters for the initial models were explored around 
the values of r c = 0.72 for the gaseous and stellar compo- 
nents and of r c = 2.2kpc with a cut off around 20kpc for 
the dark matter component. The initial inner slope of the 
density profile was chosen around 7 = 2 to produce a mas- 
sive and cuspy dark matter profile, while a more shallow 
profile with 7 = 3/2 was adopted for the gas and stellar 
profiles. 

Considered the preliminary results on the orbit anal- 
ysis from Sect. 12.41 models with an initial cuspy den- 
sity distribution for the dark matter component have 
to be preferred in order Carina to survive four pericen- 
ter passages. Models with a starting relevant amount of 
gas suffer a flattening phenomenon of the DM density 



a reduced probability to survive di fferent pericenter pas- 
sages (e.g., Penarrubia et al. ( 20101 )). 

The initial temperature for the gas can be inferred 
assuming a spherical collapse model for the initial dark 
matter model and a matter power spectrum P{k) at red- 
shift zero compatible with studies from the SDSS, e.g., 
iTegmark et al. ( 20041) or work on the analysi s of the 
Lyman-a forest, e.g.. lGnedin fc Hamilton! ( 2002 ). In this 
case the typical RMS internal velocity of a halo within 
M < 1O 8 M is less than 13km • s^ 1 and therefore we 
assumed an initial sound speed in hydrogen of roughly 
around c s = 10km- s -1 corresponding to a temperature of 
T =■ 10000°K. 

Starting with a three-component model, we implicitly 
assume our evolution to start after the oldest stellar pop- 
ulation settled to form the central nucleus of the Carina 
dwarf galaxy. Any burst of star formation that occurred 
earlier than 9-10 Gyr ago cannot help us to trace back the 
orbital parameters. Moreover, the timescale permitting to 
the initial proto-cloud to collapse, in order to form the 
initial old population, and finally to settle on the orbit, 
requires no less than 3 Gyr (according t o e.g. , the chemi- 
cal enrichment models iMarcolini et al.l (|2006|)). Thus our 
orbit determination will start on a possible Carina dwarf 
galaxy configuration of 9 Gyr ago. 



3.2. Match with the observations 

The dynamical analysis we present shows the best-fitting 
model for the data as co mpared to the observatio n of th e 
surface density profile of Irwin fc Hatzidimitriou (|l995 ). 
and more recently updated in Munoz et al.l ( 2006 ) or 
Walker et al.l ( 2007 ) with a radial velocity dispersion pro- 
file. The assumed initial and final best-fit model is pre- 
sented in Table [2] 

With the numerical scheme ad opted the number o f 
particle s is not constant (see, e.g. Pasetto et al. ( 201fj| ). 
Berczik (Il999l) . iBerczik fc Kravchukl (|200ol ) and references 
therein for detailed formulation of the star formation pro- 
cesses). The starting models are described with 350000 
particles growing with time up to half a million of parti- 
cles, and the reference model we are using for the match 
with the present day observations (after 9 Gyr of evolu- 
tion) contains about 1.2 x 10 5 stars within the 5 kpc from 
the barycentre of the Carina dwarf galaxy. 

An illustrative sketch of the evolution of the system 
after about 9.6 Gyr of evolution is presented in Fig. [7J 
The plot shows several low resolution streams where the 
mass density drops down to VlO 2 ^" VlO 3 °f ^ ne total mass 
thus showing clearly how the tidal tails developed by the 
Carina dwarf galaxy are expected to be extremely diffuse 
after 9 Gyr of evolution (i.e. now). In this sense the more 
diffuse zones have characteristics that the model can pre- 
dict but where the obser vations are lacking or t heir inter- 
pretation is insecure, e.g. iMaiewski et al.l (|2005h . The last 
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Table 2. Parameters for the galaxy model evolved in or- 
bits around the potential. Effective radius and central ve- 
locity dispersion are omitted because the whole extent of 
the profile is matched (see the text for further details) . The 
model is set at the starting time t — to in the position of 
the phase space of Eqn. (jSJ) and reaches the phase-space 
point given by Eqn. ([7]). The location in galactocentric 
coordinates is also in excellent agreement with the obser- 
vation fe.g. lMateo et all (|l998l) ) 





tlb = to 


Ub = 


1 




^ 260 [deg] 


b 




= -22 [deg] 


d 




^ 100 [kpc] 


Vl.o.s. 




219 [kms- 1 ] 


M gas 


free parameter 


2* 


M atar 


1.998 x 1O 7 M 


1.9 x 10 6 M Q 


Afdark 


1.978 x 10 8 M Q 


0.675 x 10 8 M Q 




10 



10, 



2 3 
r [kpc] 



Fig. 9. Surface density from observational data and the 
baryonic profile deduced from the observation. The model 
is in red and the data with the error bars are shown in 
blue. 
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Fig. 8. A zoom of the previous figure centred on the dwarf 
galaxy system S2. The red arrow shows the eigenvector 
direction pointing in the radial direction where the radial 
profiles are computed (see text for details). We have not 
performed our analysis in comparison with the observation 
along the l.o.s. thanks to the spherical symmetry developed 
by of the system at t=0. 



snapshots of the orbit evolution are also presented with a 
black solid line. We observe the general 'tidal tail flipping' 
phenomenon already evidence in high resolution simula- 



tion by e.g., iKlimentowski et al. (2009). The model at a 



generic snapshot is analyzed in a reference system cen- 
tered on the baryonic matter barycenter as can be seen in 
Fig.! 

We use Fig. [5J to illustrate how the tidal tails devel- 
oped in the simulation since the first pericentre passage, 
allowing us to define an ellipsoid as soon as we plot the 3D 
contours for even more diffuse regions. In Fig. [5] a yellow 
ellipsoid is visible in the center for which we can easily 
compute the inertia matrix to define its principal axis. 
Once the principal axes along the tidal tail direction are 
determined, we can perform our analysis along this radial 



direction (e.g., the red arrow in this figure). However, we 
underline that the radial stellar profiles in the orthogonal 
principal directions present similar profiles not evidencing 
any pronounced ellipsoidal shape at t = 9.6Gyr (thus pro- 
viding a good agreement w ith the observational approach 



used in lMuhoz et"~aT1 (|2006h ) 



3.2.1. Surface density profile 

Fig. P shows a plot of the density profile and the data 
ri lMun ~ 



from lMunoz et al.l (|2006|) where we have exploited the true 
physical dimensions instead of the angular dimensions for 
the system and the data representation because are an im- 
mediate indicator of the dynamical range we are consider- 
ing (e.g. the surface on the celestial sphere can be calcu- 
lated in spherical coordinates £ = dem dip J 7 dO sin 9, 
where a, /3, 7 are angular borders of the CCD mosaic of 
observational data, suitably converted from arcmin to ra- 
dian.) 

An immediately evident result of our simulations is 
the reproduction the obser yed change in the s lope of the 
surface density profile (e.g. iMunoz et al. (2006)) which fit 



the observational data quite well. 
3.2.2. Radial velocity dispersion 

The other property we can now match is the radial veloc- 
ity dispersion ay (r) which is typically considered a funda- 
mental indicator of the dynamical model of the galaxy if 
analyzed with the Jeans equation for spherical-stationary 
systems fFig. ITU]). The best- matching model nevertheless 
shows compatibility with the general trend of the baryonic 
matter traced by the observations. The general decline of 
the velocity dispersion l.bkpc is also reproduced as a fea- 
ture of the model (compare with the change of the slope of 
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Fig. 7. N-body representation of the evolved galaxy. The classical point-representation has been substituted with an 
3D isocontour plot by counting the amount of stars within a 3D grid spanning all the configuration space of the orbit 
and joining cell with equal number of stars. The green isocontours show overdensities of 1/10 4 and the purple for 
overdensities of 1/10 3 . 



the surface density profile shown in Fig. [5] The innermost 
velocity dispersion value is around 8 kms -1 . 



3.2.3. Star formation history 

At the last point we present the star formation history 
obtained from the N-body simulations. The criteria for 
star formation and the evolution of t he iso lated model 
have been adopted from lPasetto et al~ ( 2010t ). The result- 
ing SFH for the orbiting model has been normalized here 
to the the highest observational peak in the star forma- 
tion h istory deduced from the observations by iRizzi et"al 
(1200 



This approach is mainly due to the impossibility of ad- 
equately resolving a dwarf galaxy with a one-to-one cor- 
respondence between stars and particles. The true galaxy 
Carina will contain probably several million of stars and 
only a small fraction of them can be represented with an 
N-body simulation. As a consequence, =M- depends on the 
resolution of our dwarf galaxy model. The consequence of 
this limitation is that the initial amount of gas that we 



have to assign to the dwarf galaxy model for Carina when 
it starts to evolve is still a free parameter that depends on 
the gas consumption that led to the formation of its oldest 
stellar population and on the gas consumption of the self- 
consistently triggered star formation. Our goal has been 
to adapt the initial amount of gas such that there is no 
gas left after the last burst of star formation. This should 
be able to reproduce a Carina dwarf at the present time 
where essentially no gas is observed. 



While we are unable to impose strong constraints on 
this initial gas fraction and in order to reproduce the in- 
duced burst of star formation, we present in Fig. [TT] the 
simulated SFH normalized to the observed SFH. We em- 
phasize that we are not claiming to reproduce the true star 
formation history, but that we simply underline how tuned 
initial conditions are able to reproduce a system with 
present-day properties comparable with Carina including 
a star formation history dominated by well-separated star 
formation episodes. In particular, to obtain the star for- 
mation of the Carina dwarf galaxy we have to tune 
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2 3 
r [kpc] 



Fig. 10. Radial velocity dispersion profile as computed 
from the model (in red) and the observation (in blue with 
the error bars). 



120 




Fig. 11. Star Formation Rate (SFR) as function of time 
normalized to the maximum of the observed SFR (SFR Q ). 
The upper panel shows the radial distance of the N-body 
system to the origin of the reference system centered on 
the MW's barycentre for the best fitting orbit. The lower 
panel shows the star formation triggered by the pericentre 
passages and we r eported the profile of SFR as deduced in 
Rizzi et al.l (|2003). In the lower panel green and black thin 
lines refer to the gas and bound mass fraction normalized 
to the initial value (/ (t) / f in the right vertical axis). 

— a suitable amount of gas that survives the initial phase 
of SNII in the center of Carina and that initially coex- 
ists with the primordial oldest stellar population, 

— the original location of this dwarf galaxy that needs to 
be located not so far from the first pericentre passage. 
The first tidal pericentric passage may not extinguish 
the internal star formation processes. 

— the amount of dark matter at the initial condition that 
will survive the first pericentre passage stripping. 



If the simulations are started prior to the apocentre pas- 
sage (tib < 10 Gyr), the first infall into the potential well of 
MW galaxy, which takes approximately 1 Gyr, causes too 
early a gas dissipation and dark matter stripping, making 
us unable to reproduce SFH-like Carina. For this reason, 
this alternative scenario of formation of the Carina dwarf 
galaxy has not been further investigated. 

These conditions for the initial gas fraction and or- 
bital initial conditions represent our chosen initial phys- 
ical parameters to obtain the specific Carina SFH as an 
output that match the observations. As a general trend, 
we see from Fig. [TT] that the star formation after the first 
passage at the pericentre is typically more active. The 
galaxy produces an evident increase in the star formation 
rate around 6 Gyr ago and some smaller peaks correlated 
with the subsequent pericentre passages. A general trend 
is the evidence of the skewness of the star formation peaks, 
showing a rapid rise in the early phases and slower decline 
after the point of closest approach to the primary. In the 
same figure we present the evolution of the bound mass 
fraction (black line) and the gas fraction (green line). The 
initial amount of dark matter in the dwarf galaxy permits 
its baryonic component to make the first pericenter pas- 
sage almost unperturbed, shielded within the dark matter 
potential well. The dark matter result is strongly stripped 
during the first pericenter passage, while the stripping re- 
sults less important in the following pericenter passages. 
We find that gas consumption is roughly a monotonically 
decreasing function of time, controlled by the star forma- 
tion egisodes_tr^gered by the orbital pericenter passages 
fsee lRevaz et al. ( 20091 ) for a different approach based on 
completel y isolated models and different SNII recipes with 
respect to IPasetto et al. l(|2O10k 



4. Discussion and model prediction 

In this paper, we presented an extensive treatment of the 
evolution of an orbiting Carina-like dwarf galaxy and we 
self-consistently derived the star formation history. 

Efforts have been made to consider all of the available 
observational constraints. Our methodology based on the 
minimum action, and developed in order to investigate all 
possible phase-space initial conditions, leads to a family of 
orbits able to reproduce with full N-body simulation a sys- 
tem compatible with the observational constraints. Once 
the orbital parameters that minimize the action in a real- 
istic MW potential were derived, the dynamical properties 
and the star formation history of Carina were recovered. 

A previous study that attempted to interpret the ob- 
ser vational const r aints on Carina was recently presented 
byH 

unoz et al However, in this paper the au- 

thors neglected the influence of star formation processes 
in their match with the observational data. Their work 
is based on a single component model, i.e., probably, 
the dark matter, which is tuned to match the kinematic 
properties of the baryonic matter, thus fitting directly a 
one-component model to match the observed properties. 
Another important difference between our work and the 
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work by Munoz et al.l ( 2008 ) is in their use of a static ex- 
ternal potential: by increasing the size of the MW halo 
during the time range of evolution (and keeping all the 
other parameter fixed) we can vary the instants of the 
pericentre passages of Carina dwarf galaxy in dependence 
on the initial conditions allowed by the proper motion er- 
ror bars. 

In our work, the star formation history seems mainly 
driven by tidal pro cesses, contra r y to s uggestions in other 
recent works (e.g., Piatek et al. ( 20031 )). To quantify the 
influence of the environment on the star formation rate 
is always a difficult tas k not only for d warf galaxies (e.g. 



Kennicutt et"aT1 (|l994h ). Recently also iLee et all (|2009h 



showed how in the local environment a b ursty SFH ha s 
to be statisti c ally preferred (see also, e.g.. iGrebell (|19970 ; 
Mated (|l998l ); lRizzi et al.l (j2003l) h Nevertheless, we do not 
claim to have found the true reason for the episodic star 
formation history that Carina experienced in the past. We 
just presented a possible interpretation in the framework 
of the dynamical interactions. The details of the interplay 
with the environment still need to be fully understood and 
reprod uced for interacting dwarfs (e.g., Bouchard et al.l 
(|2009l) h However, by tuning the amount of gas that is col- 
lapsing in a dark matter dominated protocloud of Carina 
(and that is retained despite initial phase of Supernovae 
type II explosions), a quite similar star formation history 
for the Carina can be reproduced self-consistently as re- 
sult of a direct reaction of the system to the MW tidal 
field. 

Once other processes will emerge as clearly influencing 
the amount of gas during the history of the Carina dwarf 
galaxy, e.g., the effects of ram-pressure stripping as possi- 
ble key ingredient to clean the bu l k of i nterstellar matter 
from dwarf galaxies ([Grebel et al. I (l2003h h then the initial 
amount of gas has to be tuned again, however without in- 
fluence on our orbit determination. Nevertheless the ex- 
pected d ensity for the MW gaseous halo is probably too 
low (e.g.. lMurahl (<2000l )) in order to strongly influence the 
gas evolution of the inner part of the Carina dwarf galaxy 
which orbits between 60 and 100 kpc from the MW center 
(see, e.g., FigS]). 

Finally, the probability of an interaction with the 
Magellanic Clouds is considered in Appendix [C] 
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Appendix A: Determination of spatial velocities 

We update the tran s forma tion matr ix T d efined in 
Johnson fc Soderbloml (|l987h (following [Green] (jl985l )) to 



equinox J2000: 



-0.0534 -0.8750 -0.4810 
+0.4939 -0.4418 +0.7488 
-0.8678 -0.1975 +0.4558 



The matrix B of iJohnson fc Soderbloml (|1987l) remains 
defined as: 

cos(a) cos(<5) — sin(a) — cos(a) sm{5) 
B = T ■ | sin(o:) cos((5) cos(a) — sin(Q;) sin(<5) 
sm(5) cos((5) 

where (a,<5) are the equatorial coordinates of the generic 
dwarf galaxy at the present time t = to- In our case we get 
(a, <5) Car = (1-76, —0.89) rad for the Carina dwarf galaxy. 
Now we have to take into account a reflection of the ve- 
locity axis that we want pointing away from the Galactic 
center and we get the correction for the rotation of the 
galaxy as well as the motion of the sun relative to the 
local standard of rest: 



+ 









-1 
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V o o i) 




I V xQi LSR \ 


\\\ 


Vy0,LSR 




\ VzQ,LSR } ) ) 




where v = {v x ,v y ,v z } is the velocity of the dwarf 
galaxy at the instant t = to, v r is the radial ve- 
locity, (fi a ,fis) the proper motions in arcsec-yr -1 , 
d the distances that have to be assumed in pc to 



conversion 
l 



values 
as in 



4.74; Vc^.l sr 



Binnev fc Merrifieldl (|l998h 



apply the 
{10.0,5.2,7.2} km - s 
and V c = 220km • s" 1 as adopted by the IAU (1986). For 
the d istances, we also adopt a diff erent reference system 
from IJohnson fc Soderbloml |l987) in order to obtain a 
right-handed reference system. This is achieved by impos- 
ing a double reflection for the positive X-axis originating 
from the Galactic center and the positive Y-axis pointing 
along decreasing Galactic longitude: 



cos(&) cos(Z) 
d | cos(&) sin(Z) 
sin(fe) 



where x = (a;, y, z) is the generic position of the dwarf 
galaxy, Xq = {xq^i/q, z Q } = {8.5,0.0,0.0} is the sun's 
position assumed for simplicity to lie in the plane of the 
MW and (l,b) are the Galactic coordinates that can be 
consistently derived for a cross check from 

















-1 




-1 









v 





1 




cos(&) cos(Z) 
cos(6) sin(Z) 
sin(6) 



cos (S) cos (a) 
cos (5) sin (a) 
sin (S) 



The selection of an aligned reference system between 
the configuration space and the velocity space will permit 
us to treat, in a simpler form, the velocity vector as a 
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derivative of the position vector for times different from 
the present. Similarly, we can derive the errors as 

r 2 




P 



B22B23 
B32B33 



(I) 




k 2 + 



where p is the parallax and B 2 the matrix with elements 
{Bij) 2 With these equations and the proper motion 
and errors in the text, we can derive our estimate value for 
the orbital energy of the dwarf galaxy we are analyzing. 



Appendix B: Barycentre determination 

Once we try to move from the crude point mass determi- 
nation to the full N-body description and include many 
different astrophysical aspects, we encounter the necessity 
to analyze the system properties of the many N-body sys- 
tem realizations by detecting the position of the center of 
mass of the system. The determination of the not inertial 
center of mass of a galaxy moving within an external force 
field can be performed in different ways but it is clear that 
we must evaluate it automatically to speed up the anal- 
ysis of the large number of simulations we performed for 
finding the best match with the observational constraints. 
In this appendix, we present an original approach based 
on the downhill simplex method (Nelder & Mead 1965). 
The downhill simplex method is a direct search method 
that works moderately well in low-dimensional stochas- 
tic problems. Our task is to apply the method to find the 
barycentre of the nucleus of our dwarf galaxy, limiting our- 
selves to the knowledge of the position of the galaxy at a 
given snapshot Xi (t) £ R 3 for every particle, Vi, and at 
every moment, Vi, during the evolution of our dwarf. We 
refer the reader t o books such as Numerical Recipes by 
Press et al.l ([19861 ) which used the method in the section 
dedicated to the minimization already back in 1986. Here 
we limit ourselves to showing a further possible applica- 
tion where this method is suitable. In the practical imple- 
mentation, we refer to the matrix formalism developed in 
the original work of Nelder and Mead. The simplex, a con- 
vex hull of a tetrahedron for our 3D space surrounding our 
MW galaxy be represented by a 3x (3 + 1) time-dependent 
matrix, whose columns are the vertex 



a(t) = (B(t) x 4 (£)) , 

where B = (x\ (t) ,X2 (t) ,^3 (i)). For any simplex a £ 
M 3 we can define the matrix, M = M m as the 3x3 
matrix whose columns represent the edges of a (i) : 

M (t) = (xi — a?4, X2 — X4, X3 — X4)\ t _f — B (t) — x^e 1 " 

and e = (1, 1, 1) T . In this way we can, as a first step, check 
the degenerate character of the simplex by ensuring that 



the 3D volume of a (t), vol (a {£)) = \ |det (M (*)) | > 0. 
As a consequence, in a Euclidean geometry the reflec- 
tion, expansion, inside/outside-contraction and shrinkage 
computed by the algorithm will always produce a non- 
degenerate tetrahedron. We define then the diameter of 
the simplex as (a) = max \\xi — Xj\\ Vt where |.| is the 

standard norm. Finally we define a suitable function for 
finding the best point representing a barycentre of the 
dwarf galaxy. If we call di = ||a;M — x&|| the distance of 
the i tft, -star from the guess value of the barycentre Xb at 
the instant t = t under consideration, then we need to 
maximize the function 



di<r. 



(B.l) 



where is a characteristic radius of the system we want 
to analyze. Our experience shows that it does not have to 
be physically related to the system, but a suitable choice 
of a few kpc can much improve the convergence veloc- 
ity and the stability of the barycentre value if related to 
the convergence criteria. We know, in fact, that if / is 
a bounded function then for every non-degenerate case 
it can be proved for the downhill simplex algorithm (we 
indicate with k the iteration of the algorithm) that: 

— the sequence j/ y x i'j J always converges; 

— at every non-shrinking iterations fc, / (x[ k+1 ^ ^ 

/ f 33 ^) for 1 ^ i ^ n + 1, with strict inequality for 
at least one variable of I; 

— if there is only a finite number of shrink iterations, 
then 

— each sequence |/ (^ x ^^j j converges as k — > 00 for 
1 < i < n + 1; 

then 



— if we call 



/ [x^ ) < / (x[ k) ) Vfc and 1 < i < n + 1; 



— if there is only a finite number of non-shrink iterations, 
then all simplex vertexes converge to a single point. 

All this ensures that the down hill simplex method applied 
to the function (|B. 1[) will always converge if the number of 
shrink iteration is small (which is the case in our experi- 
ence) . We adopted a variable diameter as indicated by the 
subset of the stars for which we compute the barycentre, 



(a) < : lim 7*^ (i) = 0, e.g. with a physical radius 

3-t-oo 

converging to zero as the number of successive times that 
the down- hill simplex method is applied, j, increases the 
assignedVi = t. In practice we of course chose r* to be 
slightly higher than the softening length of the gravita- 
tional potential > e s t a r- 
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Appendix C: Encounter with the Magellanic 
Clouds 

In this appendix we discuss the probability for a fly- 
by interaction with the Ma gellanic Clouds to i nvestigate 
the possibility suggested by Muhoz et al.l ( 20061 ). We rele- 



gate these arguments to this appendix because they sim- 
ply concern probability considerations that are derived 
from dynami c s. The suggestion presented in the paper of 
iMunoz et ail (|2006l ) comes from the analysis of only 15 



stars in the local universe, too small to justify a further 
set of full simulations. In particular, we are interested in 
the recent determination of surprisingly large proper mo- 
tion for the Large Magellani c Cloud (LMC) and the Sm all 
Magellanic Cloud (SMC) (IKallivavalil et al.1 (l2006bllaR 



which caused interest due to th e possible implicat i ons fo r 
their evoluti o nary h istory Ce.g. IKallivavalil et al.l ( 2009); 



Piatek et all (l2008h: iBesla et al.l (120071); lOlsen fc Massev 



20071) ). As explained in IKallivavalil et al.l (|2008l) . the 
most relevant problem faced by the integration of the or- 
bits of the Magellanic Clouds with an expected apocentre 
so far from the MW bary centre, is the unknown mass dis- 
trib ution of the MW fo r distances larger than 100-200 kpc 
(e.g. IBesla et al.l (|2007l )). With our MW galactic potential 



model presented in Appendix B, we cannot extend the in- 
tegration of the LMC orbit much further as required to 
overlap the entire time with the orbit integration spanned 
by the minimization of the action for Carina. Thus we will 
limit ourselves to the last 3 Gyr of evolution in look-back 
time because in the potential of MW, the LMC after 3 
Gyr of look-back time integration will already be more 
distant than 200 kpc from the MW barycentre, falsifying 
any phase-space derivation. Thus we will accept the con- 



clusion already presented in IBesla et al.l (|2007l ) in favor 



of a single/first pericentre passage for the MC not more 
than few hundred million years ago and we proceed by 
minimizing the distance function dcar-MC '■ ^ 7 ^ + 



dcar-MC — dcar-MC { v 0,Car , V(),MC , t) 

in the 7-dimensional space of the initial values for Vo.car € 
[v 0<O ar ± Sv 0i car], v . M c € [v ,mc ± Svq.mc] computed 
as in Appendix A and t £ ]— 3, 0] Gyr. We integrated the 
equation of motion here for the Magellanic Clouds, taking 
into account an extra term due to the dynamical friction as 
ascribed in Eqn. © caused by the dependence of the force 
on the square of the mass of the Magellanic Clouds (e.g. 
ti lmc — 2. x 10 9 M Q ) that makes the dynamical friction 
much more relevant than in the case of Carina. The result 
is that the closest distance approach permitted between 
the two galaxies in the range of the possible phase-space 
observational errors is more than 50 kpc, completely ruling 
out any possible interaction within the currently suggested 
phase-space error range deduced from the observations. 
The same consideration holds for an interaction with the 
SMC. 



